Dublin Data Science Workshop on Bayesian Data Analysis
Mick Cooney mickcooney@gmail.com
Monday, 7 October 2019
All code and data for this workshop is available at the following URL:
https://github.com/kaybenleroll/data_workshops
Code is available in the ws_bayesian_201910/ directory.
Much of content of this workshop is based on the book “Doing Bayesian Data Analysis” by John Kruschke.
https://sites.google.com/site/doingbayesiandataanalysis/.
Another excellent book is “Data Analysis Using Regression and Multilevel/Hierarchical Models” by Gelman and Hill
http://www.stat.columbia.edu/~gelman/arm/
For a more technical and comprehensive of the area, there is “Bayesian Data Analysis” by Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin
http://www.stat.columbia.edu/~gelman/book/
Finally, I highly recommend the book ``Statistical Rethinking’’ by Richard McElreath. It teaches basic statistical principles from a Bayesian perspective.
http://xcelab.net/rm/statistical-rethinking/
1 Conditional Probability
Suppose that in the general population, the probability of having a specific rare disease (the Dreaded Lurgy) is one in a thousand. We denote the true presence or absence of the disease as the value of a parameter, \(\theta\), that takes the value 1 if disease is present, or 0 if the disease is absent. The base rate of the disease is therefore denoted \(p(\theta = 1) = 0.001\).
This is our prior belief that a person selected at random has the disease.
Suppose that there is a test for the disease that has a 99% hit rate - if a person has the disease, then the test result is positive 99% of the time.
We denote a positive test result as \(D = 1\), and a negative test result as \(D = 0\). The observed test result is a bit of data that we will use to modify our belief about the value of the underlying disease parameter.
The hit rate is expressed as
\[ p(D = 1 \, | \, \theta = 1) = 0.99. \]
The test also has a false alarm rate of 5%.
This means that 5% of the time when the disease is not present, the test falsely indicates that the disease is present. We denote the false alarm rate as
\[ p(D = 1 \, | \, \theta = 0) = 0.05 \]
However, what we need to know is \(p(\theta = 1 \, | \, D = 1)\), i.e. the probability that the patient has the disease given a positive test result.
We can calculate the above conditional probability given Bayes’ Rule and using arithmetic and algebra but we will use simulation to estimate it.
We will try to estimate this probability using simulation.
n_sim <- 1000000
base_rate <- 0.001
true_rate <- 0.99
fa_rate <- 0.05
single_test_tbl <- tibble(id = 1:n_sim) %>%
mutate(sick_person = sample(c(TRUE, FALSE),
n_sim,
prob = c(base_rate, 1 - base_rate),
replace = TRUE))
sick_tbl <- single_test_tbl %>% filter(sick_person == TRUE)
notsick_tbl <- single_test_tbl %>% filter(sick_person == FALSE)
sick_tbl <- sick_tbl %>%
mutate(test_result = sample(c(TRUE, FALSE),
n(),
prob = c(true_rate, 1 - true_rate),
replace = TRUE))
notsick_tbl <- notsick_tbl %>%
mutate(test_result = sample(c(TRUE, FALSE),
n(),
prob = c(fa_rate, 1 - fa_rate),
replace = TRUE))
single_test_tbl <- list(sick_tbl, notsick_tbl) %>%
bind_rows() %>%
arrange(id)1.1 Questions
We ask some simple questions first:
- How many people are sick?
- How many people test positive?
- Of the people with a positive result, how many are actually sick?
single_test_tbl %>% count(sick_person)## # A tibble: 2 x 2
## sick_person n
## <lgl> <int>
## 1 FALSE 999079
## 2 TRUE 921single_test_tbl %>% count(test_result)## # A tibble: 2 x 2
## test_result n
## <lgl> <int>
## 1 FALSE 948904
## 2 TRUE 51096single_test_tbl %>% count(sick_person, test_result)## # A tibble: 4 x 3
## sick_person test_result n
## <lgl> <lgl> <int>
## 1 FALSE FALSE 948896
## 2 FALSE TRUE 50183
## 3 TRUE FALSE 8
## 4 TRUE TRUE 913single_test_tbl %>%
filter(test_result == TRUE) %>%
summarise(n_sick = sum(sick_person), sick_prop = n_sick / n())## # A tibble: 1 x 2
## n_sick sick_prop
## <int> <dbl>
## 1 913 0.01791.2 Inputs Dependency
We now want to see how this probability depends on both the accuracy of the test as well as the false alarm rate.
We start looking at changing the accuracy of the test from 90% to 100%.
calc_acc_prop_func <- function(acc_rate) {
sick_prop <- create_medtest_data(100000,
base_rate = 0.001,
true_rate = acc_rate,
fa_rate = 0.05) %>%
filter(test_result == TRUE) %>%
summarise(sick_prop = sum(sick_person) / n()) %>%
pull(sick_prop)
return(sick_prop)
}
acc_props_tbl <- tibble(acc_vals = seq(0.90, 1.00, by = 0.001)) %>%
mutate(sick_prop = map_dbl(acc_vals, calc_acc_prop_func))
ggplot(acc_props_tbl) +
geom_line(aes(x = acc_vals, y = sick_prop)) +
expand_limits(y = 0) +
xlab('False Alarm Rate') +
ylab('Proportion of True Positives')
We now want to look at the effect of the false alarm rate on the conditional probability.
calc_fa_prop_func <- function(fa_rate) {
sick_prop <- create_medtest_data(100000,
base_rate = 0.001,
true_rate = 0.99,
fa_rate = fa_rate) %>%
filter(test_result == TRUE) %>%
summarise(sick_prop = sum(sick_person) / n()) %>%
pull(sick_prop)
return(sick_prop)
}
fa_props_tbl <- tibble(fa_rate = seq(0.001, 0.10, by = 0.001)) %>%
mutate(sick_prop = map_dbl(fa_rate, calc_fa_prop_func))
ggplot(fa_props_tbl) +
geom_line(aes(x = fa_rate, y = sick_prop)) +
expand_limits(y = 0) +
xlab('False Alarm Rate') +
ylab('Proportion of True Positives')
What if we have multiple tests?
fa_vals <- seq(0.001, 0.10, by = 0.001)
calc_fa_prop_func <- function(fa_rate) {
sick_prop <- create_two_medtest_data(100000,
base_rate = 0.001,
true_rate = 0.99,
fa_rate = fa_rate) %>%
filter(test_result_1 == TRUE, test_result_2 == TRUE) %>%
summarise(sick_prop = sum(sick_person) / n()) %>%
pull(sick_prop)
return(sick_prop)
}
sick_2_prop_tbl <- tibble(fa_rate = fa_vals,
sick_2_prop = map_dbl(fa_vals, calc_fa_prop_func))
sick_tbl <- fa_props_tbl %>%
inner_join(sick_2_prop_tbl, by = 'fa_rate')
ggplot(sick_tbl) +
geom_line(aes(x = fa_vals, y = sick_prop), colour = 'red') +
geom_line(aes(x = fa_vals, y = sick_2_prop)) +
expand_limits(y = 0) +
xlab('False Alarm Rate') +
ylab('Proportion of True Positives')
2 Analytical Approach
Bayesian reasoning is as old as the concept of probabilities, but has only recently started to receive a lot of attention. One likely reason for this is that, apart from a few special cases, it is mot possible to perform the calculations analytically.
In our application we have a prior distribution for our beliefs, \(p(\theta)\), and a likelihood for the data, \(p(D | \theta)\), and through use of the Chain Rule, we get the posterior distribution, \(p(\theta | D)\),
\[ p(\theta | D) = \int d\theta \, p(D | \theta) \, p(\theta). \]
In those special cases, the likelihood function has a prior and posterior distribution with the same functional form, i.e. the prior and the posterior are two members of the same `family’ of functions.
For the rest of this workshop we are going to deal with estimating the fairness of a coin, based on the result of multiple coin tosses. We define `success’ as the toss coming up Heads, and denote this probability as \(\theta\). Thus,
\[ P(y = 1 | \, \theta) = \theta \; \text{ and } \; P(y = 0 | \, \theta) = 1 - \theta. \]
We can combine the above two into a single expression:
\[ P(y | \theta) = \theta^y (1 - \theta)^{(1 - y)}. \]
2.1 The Beta Distribution
Now we consider the data \(y\) to be fixed, and consider the above as a function of \(\theta\). With this approach we call the above equation the likelihood function of \(\theta\)}.
The Bernoulli function has a conjugate prior: the Beta distribution, \(\text{Beta}(a, b)\).
R supports the beta distribution natively, via the standard grouping of functions for probability distributions: rbeta(), dbeta(), pbeta(), qbeta().
We now plot the density distribution for the Beta distribution using various combinations of parameters \(a\) and \(b\).
beta_tbl <- tribble(
~a, ~b,
1, 1,
2, 2,
10, 5,
3, 5,
)
construct_beta_tbl_func <- function(a, b) {
tibble(theta = seq(0, 1, by = 0.001)) %>%
mutate(dens = dbeta(theta, a, b))
}
beta_tbl <- beta_tbl %>%
mutate(distribution = paste0("Beta(", a, ",", b, ")")
,data = map2(a, b, construct_beta_tbl_func)
) %>%
unnest(data)
ggplot(beta_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Probability Density of Beta Distribution")
One interesting observation on the parameterisation of the Beta distribution is to observe the effect of higher magnitudes of \(a\) and \(b\) while keeping their ratio constant.
beta_tbl <- tribble(
~a, ~b,
1, 1,
2, 2,
5, 5,
10, 10,
100, 100,
)
beta_tbl <- beta_tbl %>%
mutate(distribution = paste0("Beta(", a, ",", b, ")")
,data = map2(a, b, construct_beta_tbl_func)
) %>%
unnest(data)
ggplot(beta_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Probability Density of Beta Distribution")
We see that lower values represent more uncertainty for the value of \(\theta\).
Finally, it is important to note that \(\text{Beta}(a, b)\) corresponds to \(a-1\) successes and \(b-1\) failures. To see this, we look at \(\text{Beta}(2,4)\).
beta_tbl <- construct_beta_tbl_func(2, 4) %>%
mutate(distribution = 'Beta(2,4)')
ggplot(beta_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Probability Density of Beta Distribution")
2.2 Effect of Priors
We want to investigate the influence of the prior on the posterior distribution, in particular for different amounts of data.
To start, we load data with 10 trials and 1,000 trials and use these for comparisons.
cointoss_10 <- read_rds("data/cointoss10.rds")
cointoss_1000 <- read_rds("data/cointoss1000.rds")
a_10 <- sum(cointoss_10)
b_10 <- length(cointoss_10) - a_10
a_1000 <- sum(cointoss_1000)
b_1000 <- length(cointoss_1000) - a_1000
beta_bayes_tbl <- tribble(
~a, ~b,
1, 1,
2, 2,
10, 5,
3, 5,
) %>%
mutate(distribution = paste0("Beta(", a, ",", b, ")"),
prior_data = map2(a, b, construct_beta_tbl_func),
post10_data = map2(a + a_10 , b + b_10 , construct_beta_tbl_func),
post1000_data = map2(a + a_1000, b + b_1000, construct_beta_tbl_func)
) %>%
select(-a, -b) %>%
gather('type', 'data', -distribution) %>%
unnest(data)
beta_bayes_tbl %>% glimpse()## Observations: 12,012
## Variables: 4
## $ distribution <chr> "Beta(1,1)", "Beta(1,1)", "Beta(1,1)", "Beta(1,1)", "Bet…
## $ type <chr> "prior_data", "prior_data", "prior_data", "prior_data", …
## $ theta <dbl> 0.000, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, …
## $ dens <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…We have constructed the prior and posterior distributions for different sets of parameters, so we will some plots.
2.2.1 Analysis of \(N=10\) Dataset
We start with the analysis of the smaller dataset, 10 Bernoulli trials with a success rate of 4 out of 10.
plot_tbl <- beta_bayes_tbl %>%
filter(type != 'post1000_data')
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = type)) +
facet_wrap(vars(distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Prior vs Posterior Density Comparison for 10 Coin Tosses")
We want a more direct comparison of the different density plots, so rather than plotting the priors against the corresponding posterior, we instead plot the different posterior distributions against each other.
plot_tbl <- beta_bayes_tbl %>%
filter(type == 'post10_data')
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Posterior Density Comparison Plots for N=10")
We see that the prior has a strong effect on our posterior inference as the size of the dataset is relatively small, allowing the prior to have significant influence.
2.2.2 Analysis of \(N=1,000\) Dataset
We start with the analysis of the smaller dataset, 10 Bernoulli trials with a success rate of approximately 6 out of 10.
plot_tbl <- beta_bayes_tbl %>%
filter(type != 'post10_data')
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = type)) +
facet_wrap(vars(distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Prior vs Posterior Density Comparison for 1,000 Coin Tosses")
We want a more direct comparison of the different density plots, so rather than plotting the priors against the corresponding posterior, we instead plot the different posterior distributions against each other.
plot_tbl <- beta_bayes_tbl %>%
filter(type == 'post1000_data', theta >= 0.35, theta <= 0.45)
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Posterior Density Comparison Plots for N=1,000")
3 Numerical Solutions Using a Discrete Grid
For many applications, the use of simple conjugate priors is not appropriate, and we need to deal with the posterior integral calculation itself. Since analytical solutions do not exist, we use numerical techniques to approximate the integral.
The supplied functions calc_data_prob() and calc_posterior() perform these calculations. The major benefit of this approach is that we can now use arbitrary priors and perform the integration numerically.
3.1 Posteriors for \(N=10\)
We start with the smaller dataset with \(N=10\).
3.1.1 Beta(1,1)
The \(\text{Beta}(1,1)\) uniform prior is a natural place to start and we want to calculate the numerical integral using our grid approximation to calculate the posterior density.
prior_1_1_tbl <- tibble(theta = seq(0, 1, by = 0.001)) %>%
mutate(dens = dbeta(theta, 1, 1),
dens_log = dbeta(theta, 1, 1, log = TRUE)
)
data_10_loglik_tbl <- calc_data_loglik(cointoss_10, prior_1_1_tbl$theta)
num_1_1_post_tbl <- calc_posterior(prior_1_1_tbl, data_10_loglik_tbl)
num_1_1_post_tbl %>% glimpse()## Observations: 1,001
## Variables: 9
## $ theta <dbl> 0.000, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007,…
## $ dens <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ dens_log <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ loglik <dbl> NaN, -27.6370, -24.8704, -23.2546, -22.1099, -21.2233, …
## $ unnorm_loglik <dbl> NaN, -27.6370, -24.8704, -23.2546, -22.1099, -21.2233, …
## $ width <dbl> 0.000, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,…
## $ tmp <dbl> -7.745, -7.745, -7.745, -7.745, -7.745, -7.745, -7.745,…
## $ post_loglik <dbl> NaN, -19.89202, -17.12544, -15.50960, -14.36489, -13.47…
## $ post_lik <dbl> NaN, 2.29617e-09, 3.65187e-08, 1.83767e-07, 5.77309e-07…To check that this is working correctly we now plot the numerical posterior density against those calculated by the conjugate distribution.
conjugate_tbl <- beta_bayes_tbl %>%
filter(distribution == 'Beta(1,1)', type == 'post10_data') %>%
select(theta, Conjugate = dens)
numeric_tbl <- num_1_1_post_tbl %>%
select(theta, Numeric = post_lik)
plot_tbl <- conjugate_tbl %>%
inner_join(numeric_tbl, by = 'theta') %>%
gather('label', 'value', -theta) %>%
filter(!is.na(value))
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = value, colour = label)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Posterior Distribution Beta(1,1) Comparison for Conjugate and Numeric Calculations")
3.1.2 Beta(10, 5)
The \(\text{Beta}(10,5)\) prior is alternative prior that lets us investigate priors away from a mean of \(0.5\).
prior_10_5_tbl <- tibble(theta = seq(0, 1, by = 0.001)) %>%
mutate(dens = dbeta(theta, 10, 5),
dens_log = dbeta(theta, 10, 5, log = TRUE)
)
data_10_loglik_tbl <- calc_data_loglik(cointoss_10, prior_10_5_tbl$theta)
num_10_5_post_tbl <- calc_posterior(prior_10_5_tbl, data_10_loglik_tbl)
num_10_5_post_tbl %>% glimpse()## Observations: 1,001
## Variables: 9
## $ theta <dbl> 0.000, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007,…
## $ dens <dbl> 0.00000e+00, 9.97002e-24, 5.08424e-21, 1.94673e-19, 2.5…
## $ dens_log <dbl> -Inf, -52.9625, -46.7281, -43.0830, -40.4978, -38.4936,…
## $ loglik <dbl> NaN, -27.6370, -24.8704, -23.2546, -22.1099, -21.2233, …
## $ unnorm_loglik <dbl> NaN, -80.5995, -71.5986, -66.3376, -62.6077, -59.7169, …
## $ width <dbl> 0.000, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,…
## $ tmp <dbl> -7.91681, -7.91681, -7.91681, -7.91681, -7.91681, -7.91…
## $ post_loglik <dbl> NaN, -72.6827, -63.6818, -58.4208, -54.6909, -51.8001, …
## $ post_lik <dbl> NaN, 2.71842e-32, 2.20474e-28, 4.24806e-26, 1.77025e-24…Like before, we compare the numerical approximation to the conjugate version.
conjugate_tbl <- beta_bayes_tbl %>%
filter(distribution == 'Beta(10,5)', type == 'post10_data') %>%
select(theta, Conjugate = dens)
numeric_tbl <- num_10_5_post_tbl %>%
select(theta, Numeric = post_lik)
plot_tbl <- conjugate_tbl %>%
inner_join(numeric_tbl, by = 'theta') %>%
gather('label', 'value', -theta) %>%
filter(!is.na(value))
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = value, colour = label)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Posterior Distribution Beta(10,5) Comparison for Conjugate and Numeric Calculations")
3.2 Posteriors for \(N=1,000\)
We repeat this exercise for the \(N=1,000\) dataset.
3.2.1 Beta(1,1)
prior_1_1_tbl <- tibble(theta = seq(0, 1, by = 0.001)) %>%
mutate(dens = dbeta(theta, 1, 1),
dens_log = dbeta(theta, 1, 1, log = TRUE)
)
data_1000_loglik_tbl <- calc_data_loglik(cointoss_1000, prior_1_1_tbl$theta)
num_1_1_post_tbl <- calc_posterior(prior_1_1_tbl, data_10_loglik_tbl)
num_1_1_post_tbl %>% glimpse()## Observations: 1,001
## Variables: 9
## $ theta <dbl> 0.000, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007,…
## $ dens <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ dens_log <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ loglik <dbl> NaN, -27.6370, -24.8704, -23.2546, -22.1099, -21.2233, …
## $ unnorm_loglik <dbl> NaN, -27.6370, -24.8704, -23.2546, -22.1099, -21.2233, …
## $ width <dbl> 0.000, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,…
## $ tmp <dbl> -7.745, -7.745, -7.745, -7.745, -7.745, -7.745, -7.745,…
## $ post_loglik <dbl> NaN, -19.89202, -17.12544, -15.50960, -14.36489, -13.47…
## $ post_lik <dbl> NaN, 2.29617e-09, 3.65187e-08, 1.83767e-07, 5.77309e-07…conjugate_tbl <- beta_bayes_tbl %>%
filter(distribution == 'Beta(1,1)', type == 'post1000_data') %>%
select(theta, Conjugate = dens)
numeric_tbl <- num_1_1_post_tbl %>%
select(theta, Numeric = post_lik)
plot_tbl <- conjugate_tbl %>%
inner_join(numeric_tbl, by = 'theta') %>%
gather('label', 'value', -theta) %>%
filter(!is.na(value))
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = value, colour = label)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Posterior Distribution Beta(1,1) Comparison for Conjugate and Numeric Calculations")
3.2.2 Beta(10, 5)
prior_10_5_tbl <- tibble(theta = seq(0, 1, by = 0.001)) %>%
mutate(dens = dbeta(theta, 10, 5),
dens_log = dbeta(theta, 10, 5, log = TRUE)
)
data_1000_loglik_tbl <- calc_data_loglik(cointoss_1000, prior_10_5_tbl$theta)
num_10_5_post_tbl <- calc_posterior(prior_10_5_tbl, data_1000_loglik_tbl)
num_10_5_post_tbl %>% glimpse()## Observations: 1,001
## Variables: 9
## $ theta <dbl> 0.000, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007,…
## $ dens <dbl> 0.00000e+00, 9.97002e-24, 5.08424e-21, 1.94673e-19, 2.5…
## $ dens_log <dbl> -Inf, -52.9625, -46.7281, -43.0830, -40.4978, -38.4936,…
## $ loglik <dbl> NaN, -2763.70, -2487.04, -2325.46, -2210.99, -2122.33, …
## $ unnorm_loglik <dbl> NaN, -2816.66, -2533.77, -2368.54, -2251.49, -2160.83, …
## $ width <dbl> 0.000, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,…
## $ tmp <dbl> -677.315, -677.315, -677.315, -677.315, -677.315, -677.…
## $ post_loglik <dbl> NaN, -2139.350, -1856.458, -1691.228, -1574.172, -1483.…
## $ post_lik <dbl> NaN, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…conjugate_tbl <- beta_bayes_tbl %>%
filter(distribution == 'Beta(10,5)', type == 'post1000_data') %>%
select(theta, Conjugate = dens)
numeric_tbl <- num_10_5_post_tbl %>%
select(theta, Numeric = post_lik)
plot_tbl <- conjugate_tbl %>%
inner_join(numeric_tbl, by = 'theta') %>%
gather('label', 'value', -theta) %>%
filter(!is.na(value))
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = value, colour = label)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Posterior Distribution Beta(10,5) Comparison for Conjugate and Numeric Calculations")
3.3 Calculating Posteriors for More Complex Priors
Suppose we know that the trial is biased 3/1, but we do not know the direction of the bias.
To construct this prior, we average out the both Beta distributions as our prior and then feed this into out posterior calculation.
biased_prior_tbl <- tibble(
type = 'prior_data',
theta = seq(0, 1, by = 0.001)
)
beta_low_tbl <- biased_prior_tbl %>%
mutate(distribution = 'Beta (2,4)',
dens = dbeta(theta, 2, 4))
beta_upr_tbl <- biased_prior_tbl %>%
mutate(distribution = 'Beta (4,2)',
dens = dbeta(theta, 4, 2))
biased_prior_tbl <- list(beta_low_tbl, beta_upr_tbl) %>%
bind_rows() %>%
group_by(theta) %>%
summarise(dens = mean(dens)) %>%
mutate(distribution = 'Compound Prior',
type = 'prior_data',
dens_log = log(dens))
plot_tbl <- list(
beta_low_tbl,
beta_upr_tbl,
biased_prior_tbl
) %>%
bind_rows()
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Density Plot for Biased Prior Distribution")
Now that we have our prior, we plug in our integration calculator and investigate the calculated posterior.
3.3.1 10-Toss Data
For the cointoss_10 data, we use our likelihood values to calculate the posterior.
data_10_loglik_tbl <- calc_data_loglik(cointoss_10, prior_1_1_tbl$theta)
biased_10data_tbl <- calc_posterior(biased_prior_tbl, data_10_loglik_tbl) %>%
transmute(theta, type = 'posterior', post_loglik, post_lik)
plot_tbl <- list(
biased_prior_tbl %>% transmute(theta, type, dens),
biased_10data_tbl %>% transmute(theta, type, dens = post_lik),
beta_bayes_tbl %>%
filter(distribution == 'Beta(2,2)', type == 'post10_data') %>%
transmute(theta, type = 'Beta(2,2) posterior', dens)
) %>%
bind_rows()
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = type)) +
xlab(expression(theta)) +
ylab("Posterior Density")
3.3.2 1000-Toss Data
We also perform this comparison for 1,000 data points.
data_1000_loglik_tbl <- calc_data_loglik(cointoss_1000, prior_1_1_tbl$theta)
biased_1000data_tbl <- calc_posterior(biased_prior_tbl, data_1000_loglik_tbl) %>%
transmute(theta, type = 'posterior', post_loglik, post_lik)
plot_tbl <- list(
biased_prior_tbl %>% transmute(theta, type, dens),
biased_1000data_tbl %>% transmute(theta, type, dens = post_lik),
beta_bayes_tbl %>%
filter(distribution == 'Beta(2,2)', type == 'post1000_data') %>%
transmute(theta, type = 'Beta(2,2) posterior', dens)
) %>%
bind_rows()
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = type)) +
xlab(expression(theta)) +
ylab("Posterior Density")
3.4 Strong Priors
The previous prior ended up concentrating probability mass in the middle of the distribution.
We repeat this experiment now with strong priors - first observing the prior itself.
biased_prior_tbl <- tibble(
type = 'prior_data',
theta = seq(0, 1, by = 0.001)
)
strong_prior_tbl <- construct_additive_prior(biased_prior_tbl, c(20, 40), c(40, 20))
ggplot(strong_prior_tbl) +
geom_line(aes(x = theta, y = dens, colour = distribution)) +
xlab(expression(theta)) +
ylab("Density") +
ggtitle("Density Plot for Biased Prior Distribution")
We now use our numerical integration routines to calculate the posterior distribution due to this prior on the observation of 10 datapoints.
use_prior_tbl <- strong_prior_tbl %>%
filter(distribution == 'Compound Prior')
data_10_loglik_tbl <- calc_data_loglik(cointoss_10, prior_1_1_tbl$theta)
strong_10data_tbl <- calc_posterior(use_prior_tbl, data_10_loglik_tbl) %>%
transmute(theta, type = 'posterior', post_loglik, post_lik)
plot_tbl <- list(
use_prior_tbl %>%
transmute(theta, type, dens),
strong_10data_tbl %>%
transmute(theta, type, dens = post_lik),
beta_bayes_tbl %>%
filter(distribution == 'Beta(2,2)', type == 'post10_data') %>%
transmute(theta, type = 'Beta(2,2) posterior', dens)
) %>%
bind_rows()
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = type)) +
xlab(expression(theta)) +
ylab("Posterior Density")
We also contrast this with the 1,000 datapoint data.
use_prior_tbl <- strong_prior_tbl %>%
filter(distribution == 'Compound Prior')
data_1000_loglik_tbl <- calc_data_loglik(cointoss_1000, prior_1_1_tbl$theta)
strong_1000data_tbl <- calc_posterior(use_prior_tbl, data_1000_loglik_tbl) %>%
transmute(theta, type = 'posterior', post_loglik, post_lik)
plot_tbl <- list(
use_prior_tbl %>%
transmute(theta, type, dens),
strong_1000data_tbl %>%
transmute(theta, type, dens = post_lik),
beta_bayes_tbl %>%
filter(distribution == 'Beta(2,2)', type == 'post1000_data') %>%
transmute(theta, type = 'Beta(2,2) posterior', dens)
) %>%
bind_rows()
ggplot(plot_tbl) +
geom_line(aes(x = theta, y = dens, colour = type)) +
xlab(expression(theta)) +
ylab("Posterior Density")
4 First Stan Program
We now turn our attention to writing Stan programs, and showing how we can use this flexible probabilistic programming language to model all sorts of complex probabilistic models.
4.1 Parameterising the Beta Distribution
We have discussed how the Beta distribution, \(\text{Beta}(a, b)\), with two shape parameters \(a\) and \(b\) - \(a\) corresponding to the count of successful trials, and \(b\) the failed trials.
While a useful parameterisation, we can convert these variables to \((\mu, K)\).
In this model, we have
\[\begin{eqnarray*} a &=& \mu K \\ b &=& (1 - \mu) K \end{eqnarray*}\]
In this parameterisation, \(\mu\) is the underlying mean of the distribution and \(K\) is a measure of the precision of the distribution around this value - the higher the value of \(K\), the tighter the distribution is around \(\mu\).
4.2 Binomial Likelihood with Beta Prior
We now look at some Stan code to sample the posterior of a model with a Beta prior and a binomial likelihood.
data {
real prior_sh1;
real prior_sh2;
real K;
int N;
int y[N];
int prior_pd;
}
parameters {
real<lower=0,upper=1> theta;
real<lower=0,upper=1> mu;
}
model {
real a;
real b;
a = mu * K;
b = (1 - mu) * K;
mu ~ beta(prior_sh1, prior_sh2);
theta ~ beta(a, b);
if(prior_pd == 1)
y ~ bernoulli(theta);
}
generated quantities {
int y_check;
y_check = bernoulli_rng(theta);
}To use Stan, we need to build the model - the Stan code is compiled into a C++ program and the model is conditioned on the data we provide.
stan_seed <- 42
standata_lst <- list(
prior_sh1 = 2,
prior_sh2 = 2,
K = 5,
N = length(cointoss_10),
y = cointoss_10
)
model_01_stanmodel <- stan_model("first_stan_model.stan")
model_01_stanfit <- sampling(
object = model_01_stanmodel,
data = standata_lst %>% list_modify(prior_pd = 1),
iter = 1000,
chains = 4,
seed = stan_seed
)
plot(model_01_stanfit, pars = c('theta', 'mu')) +
xlab("Value") +
ylab("Parameter") +
ggtitle("Posterior Plots of the Model Parameters")
5 Checking Fit Validity
We need to make sure that the posterior we explore and sample with Stan is valid - improper specification of the model will result in invalid inferences.
This problem is real, and needs to be checked against. Unfortunately there are no methods to prove a given sample is valid, so we instead rely on heuristics and diagnostics.
A first check is to look for divergences, tree-depth issues and the energy diagnostic. We want no problems with this as a starting point.
model_01_stanfit %>% check_hmc_diagnostics()##
## Divergences:## 0 of 2000 iterations ended with a divergence.##
## Tree depth:## 0 of 2000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.Another useful diagnostic plot is to visually inspect the traceplots - a line plot of each parameter value against the index number of each iteration.
A valid sample means each parameter looks like a ‘hairy caterpillar’.
traceplot(model_01_stanfit, pars = c('mu', 'theta'), nrow = 2) +
ggtitle("Traceplots of Model Parameters")
5.1 Predictive Model Checks
A key component to checking model validity is through the use of predictive checks. There are two main types of predictive check in this context:
- prior predictive checks, where we do not use the data to condition the model - allowing us to check that our probabilistic model is sensible regardless of the data we observed.
- posterior predictive checks, where we condition the model with the data and make sure our model outputs are consistent with what we observed from the data.
Posterior predictive checks are more common, but both are useful for model checking.
We use prior checks to ensure that our models are reasonable and produce values that are valid - we think of it as illustrating the possible outputs of the model built, regardless of the data observed.
model_01_prior_stanfit <- sampling(
object = model_01_stanmodel,
data = standata_lst %>% list_modify(prior_pd = 0),
iter = 1000,
chains = 4,
seed = stan_seed
)
model_01_prior_stanfit %>% plot(pars = c('theta', 'mu'))
We use the function tidy_draws() to combine the outputs of the modelling and compare the plots of the prior and posterior parameter distributions.
prior_post_bern_tbl <- list(
prior = model_01_prior_stanfit %>% tidy_draws(),
posterior = model_01_stanfit %>% tidy_draws()
) %>%
bind_rows(.id = 'distribution')
prior_post_bern_tbl %>% glimpse()## Observations: 4,000
## Variables: 14
## $ distribution <chr> "prior", "prior", "prior", "prior", "prior", "prior", "…
## $ .chain <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ .iteration <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ .draw <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ theta <dbl> 0.540693, 0.604171, 0.706733, 0.579115, 0.187584, 0.292…
## $ mu <dbl> 0.747482, 0.440184, 0.365068, 0.332971, 0.337234, 0.537…
## $ y_check <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1…
## $ lp__ <dbl> -4.80798, -3.93989, -5.07020, -4.54789, -4.23168, -4.29…
## $ accept_stat__ <dbl> 0.9719015, 0.9669799, 0.6454783, 0.9786070, 0.9702054, …
## $ stepsize__ <dbl> 0.520804, 0.520804, 0.520804, 0.520804, 0.520804, 0.520…
## $ treedepth__ <dbl> 2, 2, 1, 1, 2, 1, 3, 2, 3, 2, 3, 1, 1, 2, 2, 2, 3, 2, 2…
## $ n_leapfrog__ <dbl> 7, 7, 1, 3, 7, 3, 15, 3, 7, 5, 7, 3, 1, 3, 3, 7, 7, 3, …
## $ divergent__ <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ energy__ <dbl> 8.95569, 5.07271, 5.14263, 5.27711, 5.27957, 5.64189, 7…We produce a number of comparison plots for \(\mu\) first:
ggplot(prior_post_bern_tbl) +
geom_histogram(aes(x = mu, fill = distribution),
bins = 30, alpha = 0.5, position = 'identity') +
xlab(expression(mu)) +
ylab("Frequency") +
ggtitle(expression(paste("Comparison Histograms for Prior and Posterior ",
"Distributions for ", mu)))
plot_tbl <- prior_post_bern_tbl %>%
group_by(distribution) %>%
summarise(low_80 = quantile(mu, 0.10),
low_50 = quantile(mu, 0.25),
median = median(mu),
mean = mean(mu),
upr_50 = quantile(mu, 0.75),
upr_80 = quantile(mu, 0.90)
)
ggplot(plot_tbl, aes(x = distribution)) +
geom_errorbar(aes(ymin = low_80, ymax = upr_80), width = 0, size = 1) +
geom_errorbar(aes(ymin = low_50, ymax = upr_50), colour = 'red', width = 0, size = 2) +
geom_point(aes(y = mean), colour = 'blue', size = 3) +
geom_point(aes(y = median), size = 3) +
expand_limits(y = 0) +
xlab("Distribution") +
ylab(expression(mu)) +
ggtitle(expression(paste("Comparison 50% and 80% Intervals for Prior and ",
"Posterior Comparisons of ", mu)))
And now we look at similar plots for \(\theta\).
ggplot(prior_post_bern_tbl) +
geom_histogram(aes(x = theta, fill = distribution),
bins = 30, alpha = 0.5, position = 'identity') +
xlab(expression(theta)) +
ylab("Frequency") +
ggtitle(expression(paste("Comparison Histograms for Prior and Posterior ",
"Distributions for ", theta)))
plot_tbl <- prior_post_bern_tbl %>%
group_by(distribution) %>%
summarise(low_80 = quantile(theta, 0.10),
low_50 = quantile(theta, 0.25),
median = median(theta),
mean = mean(theta),
upr_50 = quantile(theta, 0.75),
upr_80 = quantile(theta, 0.90)
)
ggplot(plot_tbl, aes(x = distribution)) +
geom_errorbar(aes(ymin = low_80, ymax = upr_80), width = 0, size = 1) +
geom_errorbar(aes(ymin = low_50, ymax = upr_50), colour = 'red', width = 0, size = 2) +
geom_point(aes(y = mean), colour = 'blue', size = 3) +
geom_point(aes(y = median), size = 3) +
expand_limits(y = 0) +
xlab("Distribution") +
ylab(expression(theta)) +
ggtitle(expression(paste("Comparison 50% and 80% Intervals for Prior and ",
"Posterior Comparisons of ", theta)))
5.2 Remodelling with Binomial Likelihood
While using Bernoulli distributions to model binary outcomes, it is better to compress our data into counts of successes rather than enumerating them separately.
This results in much faster computation so we can build bigger and more sophisticated models if we need to.
t10_count <- sum(cointoss_10)
t1000_count <- sum(cointoss_1000)We now look at our Stan Binomial model.
data {
real prior_sh1;
real prior_sh2;
real K;
int n;
int k;
int prior_pd;
}
parameters {
real<lower=0,upper=1> theta;
real<lower=0,upper=1> mu;
}
model {
real a;
real b;
a = mu * K;
b = (1 - mu) * K;
mu ~ beta(prior_sh1, prior_sh2);
theta ~ beta(a, b);
if(prior_pd == 1)
k ~ binomial(n, theta);
}
generated quantities {
int k_ppd;
k_ppd = binomial_rng(n, theta);
}We now construct a dataset to pass into our Binomial model and condition our model on it.
The results of this model should be the same as from our Bernoulli version (allowing for small variances due to noise).
standata_binom_lst <- list(
prior_sh1 = 2,
prior_sh2 = 2,
K = 5,
n = 10,
k = t10_count
)
model_02_stanmodel <- stan_model("simple_binomial.stan")
model_02_stanfit <- sampling(
object = model_02_stanmodel,
data = standata_binom_lst %>% list_modify(prior_pd = 1),
iter = 1000,
chains = 4,
seed = stan_seed
)
model_02_prior_stanfit <- sampling(
object = model_02_stanmodel,
data = standata_binom_lst %>% list_modify(prior_pd = 0),
iter = 1000,
chains = 4,
seed = stan_seed
)
model_02_stanfit %>% check_hmc_diagnostics()##
## Divergences:
##
## Tree depth:
##
## Energy:traceplot(model_02_stanfit, pars = c('mu', 'theta'), nrow = 2) +
ggtitle("Traceplots of Binomial Model Parameters")
plot(model_02_stanfit, pars = c('theta', 'mu')) +
xlab("Value") +
ylab("Parameter") +
ggtitle("Posterior Plots of the Binomial Model Parameters")
prior_post_binom_tbl <- list(
prior = model_02_prior_stanfit %>% tidy_draws(),
posterior = model_02_stanfit %>% tidy_draws()
) %>%
bind_rows(.id = 'distribution')
prior_post_bern_tbl %>% glimpse()## Observations: 4,000
## Variables: 14
## $ distribution <chr> "prior", "prior", "prior", "prior", "prior", "prior", "…
## $ .chain <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ .iteration <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ .draw <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ theta <dbl> 0.540693, 0.604171, 0.706733, 0.579115, 0.187584, 0.292…
## $ mu <dbl> 0.747482, 0.440184, 0.365068, 0.332971, 0.337234, 0.537…
## $ y_check <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1…
## $ lp__ <dbl> -4.80798, -3.93989, -5.07020, -4.54789, -4.23168, -4.29…
## $ accept_stat__ <dbl> 0.9719015, 0.9669799, 0.6454783, 0.9786070, 0.9702054, …
## $ stepsize__ <dbl> 0.520804, 0.520804, 0.520804, 0.520804, 0.520804, 0.520…
## $ treedepth__ <dbl> 2, 2, 1, 1, 2, 1, 3, 2, 3, 2, 3, 1, 1, 2, 2, 2, 3, 2, 2…
## $ n_leapfrog__ <dbl> 7, 7, 1, 3, 7, 3, 15, 3, 7, 5, 7, 3, 1, 3, 3, 7, 7, 3, …
## $ divergent__ <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ energy__ <dbl> 8.95569, 5.07271, 5.14263, 5.27711, 5.27957, 5.64189, 7…distrib_compare_tbl <- list(
Bernouilli = prior_post_bern_tbl,
Binomial = prior_post_binom_tbl
) %>%
bind_rows(.id = 'likelihood')
distrib_compare_tbl %>% glimpse()## Observations: 8,000
## Variables: 16
## $ likelihood <chr> "Bernouilli", "Bernouilli", "Bernouilli", "Bernouilli",…
## $ distribution <chr> "prior", "prior", "prior", "prior", "prior", "prior", "…
## $ .chain <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ .iteration <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ .draw <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ theta <dbl> 0.540693, 0.604171, 0.706733, 0.579115, 0.187584, 0.292…
## $ mu <dbl> 0.747482, 0.440184, 0.365068, 0.332971, 0.337234, 0.537…
## $ y_check <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1…
## $ lp__ <dbl> -4.80798, -3.93989, -5.07020, -4.54789, -4.23168, -4.29…
## $ accept_stat__ <dbl> 0.9719015, 0.9669799, 0.6454783, 0.9786070, 0.9702054, …
## $ stepsize__ <dbl> 0.520804, 0.520804, 0.520804, 0.520804, 0.520804, 0.520…
## $ treedepth__ <dbl> 2, 2, 1, 1, 2, 1, 3, 2, 3, 2, 3, 1, 1, 2, 2, 2, 3, 2, 2…
## $ n_leapfrog__ <dbl> 7, 7, 1, 3, 7, 3, 15, 3, 7, 5, 7, 3, 1, 3, 3, 7, 7, 3, …
## $ divergent__ <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ energy__ <dbl> 8.95569, 5.07271, 5.14263, 5.27711, 5.27957, 5.64189, 7…
## $ k_ppd <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…ggplot(distrib_compare_tbl) +
geom_boxplot(aes(x = distribution, y = theta, fill = likelihood),
position = position_dodge(width = 0.5), width = 0.3) +
xlab("Distribution Type") +
ylab(expression(theta)) +
ggtitle("Comparison Boxplot of Prior and Posterior Distributions",
subtitle = 'Bernouilli and Binomial Likelihoods')
6 Modelling Count Data
So far we have focused on the simplest of probability models - binary outcome data where we need to estimate a proportion or probability.
We now extend this concept to count data - we estimate a frequency rate of some event per given unit.
7 R Environment
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 3.6.1 (2019-07-05)
## os Ubuntu 19.04
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate C.UTF-8
## ctype C.UTF-8
## tz Europe/Dublin
## date 2019-11-24
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## arrayhelpers 1.0-20160527 2016-05-28 [1] CRAN (R 3.6.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0)
## backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.1)
## bayesplot * 1.7.0 2019-05-23 [1] CRAN (R 3.6.1)
## broom 0.5.2 2019-04-07 [1] CRAN (R 3.6.0)
## callr 3.3.2 2019-09-22 [1] CRAN (R 3.6.1)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.6.0)
## cli 1.1.0 2019-03-19 [1] CRAN (R 3.6.0)
## coda 0.19-3 2019-07-05 [1] CRAN (R 3.6.1)
## codetools 0.2-16 2018-12-24 [3] CRAN (R 3.6.0)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0)
## conflicted * 1.0.4 2019-06-21 [1] CRAN (R 3.6.1)
## cowplot * 1.0.0 2019-07-11 [1] CRAN (R 3.6.1)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0)
## DBI 1.0.0 2018-05-02 [1] CRAN (R 3.6.0)
## dbplyr 1.4.2 2019-06-17 [1] CRAN (R 3.6.0)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.0)
## devtools 2.2.1 2019-09-24 [1] CRAN (R 3.6.1)
## digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.1)
## dplyr * 0.8.3 2019-07-04 [1] CRAN (R 3.6.1)
## ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.1)
## evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0)
## fansi 0.4.0 2018-10-05 [1] CRAN (R 3.6.0)
## farver 2.0.1 2019-11-13 [1] CRAN (R 3.6.1)
## forcats * 0.4.0 2019-02-17 [1] CRAN (R 3.6.0)
## fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0)
## generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.0)
## ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.6.1)
## ggridges 0.5.1 2018-09-27 [1] CRAN (R 3.6.0)
## ggstance 0.3.3 2019-08-19 [1] CRAN (R 3.6.1)
## glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 3.6.0)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0)
## haven 2.2.0 2019-11-08 [1] CRAN (R 3.6.1)
## hms 0.5.2 2019-10-30 [1] CRAN (R 3.6.1)
## htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.1)
## httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.1)
## inline 0.3.15 2018-05-18 [1] CRAN (R 3.6.0)
## jsonlite 1.6 2018-12-07 [1] CRAN (R 3.6.0)
## knitr 1.26 2019-11-12 [1] CRAN (R 3.6.1)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.6.0)
## lattice 0.20-38 2018-11-04 [3] CRAN (R 3.6.0)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0)
## lifecycle 0.1.0 2019-08-01 [1] CRAN (R 3.6.1)
## loo 2.1.0 2019-03-13 [1] CRAN (R 3.6.1)
## lubridate 1.7.4 2018-04-11 [1] CRAN (R 3.6.0)
## magrittr * 1.5 2014-11-22 [1] CRAN (R 3.6.0)
## matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.1)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.0)
## modelr 0.1.5 2019-08-08 [1] CRAN (R 3.6.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0)
## nlme 3.1-141 2019-08-01 [3] CRAN (R 3.6.1)
## pillar 1.4.2 2019-06-29 [1] CRAN (R 3.6.1)
## pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.1)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.1)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.0)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.6.0)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.6.0)
## processx 3.4.1 2019-07-18 [1] CRAN (R 3.6.1)
## ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.0)
## purrr * 0.3.3 2019-10-18 [1] CRAN (R 3.6.1)
## R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.1)
## Rcpp 1.0.3 2019-11-08 [1] CRAN (R 3.6.1)
## readr * 1.3.1 2018-12-21 [1] CRAN (R 3.6.0)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 3.6.0)
## remotes 2.1.0 2019-06-24 [1] CRAN (R 3.6.1)
## reprex 0.3.0 2019-05-16 [1] CRAN (R 3.6.0)
## rlang 0.4.2 2019-11-23 [1] CRAN (R 3.6.1)
## rmarkdown 1.17 2019-11-13 [1] CRAN (R 3.6.1)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.0)
## rstan * 2.19.2 2019-07-09 [1] CRAN (R 3.6.1)
## rstudioapi 0.10 2019-03-19 [1] CRAN (R 3.6.0)
## rvest 0.3.5 2019-11-08 [1] CRAN (R 3.6.1)
## scales * 1.1.0 2019-11-18 [1] CRAN (R 3.6.1)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0)
## StanHeaders * 2.19.0 2019-09-07 [1] CRAN (R 3.6.1)
## stringi 1.4.3 2019-03-12 [1] CRAN (R 3.6.0)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.6.0)
## svUnit 0.7-12 2014-03-02 [1] CRAN (R 3.6.0)
## testthat 2.3.0 2019-11-05 [1] CRAN (R 3.6.1)
## tibble * 2.1.3 2019-06-06 [1] CRAN (R 3.6.0)
## tidybayes * 1.1.0 2019-06-02 [1] CRAN (R 3.6.1)
## tidyr * 1.0.0 2019-09-11 [1] CRAN (R 3.6.1)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.6.0)
## tidyverse * 1.3.0 2019-11-21 [1] CRAN (R 3.6.1)
## usethis 1.5.1 2019-07-04 [1] CRAN (R 3.6.1)
## utf8 1.1.4 2018-05-24 [1] CRAN (R 3.6.0)
## vctrs 0.2.0.9003 2019-10-08 [1] Github (r-lib/vctrs@bc20422)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0)
## xfun 0.11 2019-11-12 [1] CRAN (R 3.6.1)
## xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.1)
## yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0)
## zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.6.0)
##
## [1] /usr/local/lib/R/site-library
## [2] /usr/lib/R/site-library
## [3] /usr/lib/R/library